(library (grid-gen phys)

(export CGS CGS-lam
        2CGS 2CGS-lam
        MATH:Pi

        PHYS:gam
        PHYS:Rgas
        PHYS:G
        
        ASTRO:Rsun
        ASTRO:Msun
        ASTRO:OmegaK
        ASTRO:VK
        
        set-gam
        set-scal
)

(import (guile)
        (rnrs)
        (grid-gen main)
        (grid-gen vars)
        (grid-gen code)
        (grid-gen math)
        (grid-gen symath)
)

(define g-scal (comp-float "g_scal" 1.0 #t))
(define cm-scal (comp-float "cm_scal" 1.0 #t))
(define s-scal (comp-float "s_scal" 1.0 #t))
(define K-scal (comp-float "K_scal" 1.0 #t))

(define units
 `((g   ,g-scal)
   (cm  ,cm-scal)
   (s   ,s-scal)
   (K   ,K-scal)
   
   (dyn (/ (* g cm ) (** s 2.0)))
   (erg (/ (* g (** cm 2.0)) (** s 2.0)))
   (rho (/ g (** cm 3.0)))
   (mol g)
  )
)

(define units-h (assoc->hash units))

(define scal-h (make-hash-table))

(define set-scal (lambda (u s)
  (hash-set! scal-h u s)
  ;(let ((v (hash-ref units-h u #f))
  ;      (s1 (simp-ex (exp-v-def s #t)))
  ;     )
  ;  (if (is-var? v)
  ;      (begin
  ;        (format #t "S: ~a\n" s1)
  ;        (v 'set-def! s1)
  ;       ;; (v 'set-def-expand! #f)
  ;      )
  ;      (int-err "Invalid unit in set-scal!")
  ;  )
  ;)
))

(define def-scals
  `((,g-scal 1.0)
    (,cm-scal 1.0)
    (,s-scal 1.0)
    (,K-scal 1.0)
   )
)

(define make-scals (lambda (cod)
  (let* ((uscls (hash-map->list 
                  (lambda (u s)
                    `(,(u2v (hash-ref units-h u #f)) 
                      ,(simp-ex (exp-v-def s #f))
                     )
                  )
                  scal-h
                )
         )
         (dscls (filter (lambda (sd)
                          (not (hash-ref scal-h (car sd) #f))
                        )
                        def-scals
                )
         )
         (lu (length uscls))
         (ld (length def-scals))
         (check (if (> lu ld) (int-err "Too many scals!")))
         (scls (append uscls (list-head (reverse dscls) (- ld lu))))
;;         (dbg (for-each (lambda (x) (format #t "EQ: ~a = ~a\n" (infix-out (car x)) (infix-out (cadr x)))) scls))
         (scv (apply math-cvector
                     (map (lambda (su) `(- ,(car su) ,(cadr su))) scls)
              )
         )
         (vrs (apply math-cvector (map car def-scals)))
;;         (dbg1 (map (lambda (v) (format #t "VR: ~a\n" (infix-out v))) (arr-elements vrs)))
         (slv (solve-sys scv vrs))
        )
    (hash-for-each
      (lambda (vr eq)
        (if (is-var? vr)
            (begin
              ;; (format #t "S: ~a = ~a\n" (infix-out vr) (infix-out (simp-ex eq)))
              (vr 'set-def! (simp-ex (exp-v-def (norm-form eq) #t)))
              ;; (vr 'set-def-expand! #f)
            )
            (int-err "not a variable: ~a\n" vr)
        )
      )
      slv
    )
  )
  cod
))

(define u2v (lambda (ex)
  (cond ((number? ex) ex)
        ((is-var? ex) ex)
        ((symbol? ex)
         (let ((s (hash-ref units-h ex #f)))
           (if s (u2v s) (int-err "Invalid unit in CGS ex: ~a\n" ex))
         )
        )
        ((and (list? ex) (not (null? ex)))
         (let ((op (car ex))
               (args (map u2v (cdr ex)))
              )
           (cons op args)
         )
        )
        (#t (int-err "Invalid CGS expr: ~a\n" ex))
  )
))

(define CGS (lambda (v u)
  `(/ ,v ,(u2v u))
))

(define 2CGS (lambda (v u)
  `(* ,v ,(u2v u))
))

(define 2CGS-lam (lambda (u)
  (lambda args (map (lambda (x) (2CGS x u)) args))
))
(define CGS-lam (lambda (u)
  (lambda args (map (lambda (x) (CGS x u)) args))
))

(define set-gam (lambda (v) (set! PHYS:gam (comp-float "gam" v))))

(define MATH:Pi (comp-float "pi" 3.1415926))

(define PHYS:gam  (comp-float "gam"  `(/ 5.0 3.0)))
(define PHYS:Rgas (comp-float "rgas"  (CGS 8.3143e7 '(/ erg (* K mol)))))
(define PHYS:G    (comp-float "G" 
                              (CGS 6.67e-8 '(/ (** cm 3.0) (* g (** s 2.0))))
                  )
)

(define ASTRO:Rsun (comp-float "Rsun" (CGS 6.96e10 'cm)))
(define ASTRO:Msun (comp-float "Msun" (CGS 2.0e33 'g)))
(define ASTRO:OmegaK (lambda (m a) `(sqrt (/ (* ,PHYS:G ,m) (** ,a 3.0)))))
(define ASTRO:VK (lambda (m a) `(sqrt (/ (* ,PHYS:G ,m) ,a))))

)

(add-wrapper 7 'RUN_CODE make-scals)
